#include <Configuration.h>
#include <Utilities/RandomGenerator.h>
#include <mpi/communicator.h>
#include <mpi/collectives.h>
#include <Message/OpenMP.h>
#include <Utilities/Timer.h>
#include <OhmmsPETE/OhmmsVector.h>
#include <benchmark/vorticity.h>

/** vorticity operations
 *
 */
int main(int argc, char** argv)
{
  using namespace APPNAMESPACE;
  using namespace std;
  mpi::environment env(argc,argv);
  mpi::communicator mycomm;
  OhmmsInfo ohmms("vor",mycomm.rank(),0,1);

  typedef double real_type;
  typedef std::complex<double> complex_type;
  typedef Vector<real_type>  container_type;

  const int n_var=3;
  int niters=10;
  int n=12288;
  int npx=512;
  int npy=512;
  bool debug=false;

  int ic=0;
  while(ic<argc)
  {
    string a(argv[ic]);
    if(a.find("-d")<a.size()) //get the z
      n=atoi(argv[++ic]);
    if(a.find("-npx")<a.size()) //x-partition
      npx=atoi(argv[++ic]);
    if(a.find("-npy")<a.size()) //x-partition
      npy=atoi(argv[++ic]);
    else if(a.find("-i")<a.size()) 
      niters=atoi(argv[++ic]);
    else if(a.find("-debug")<a.size())  
      debug=true;
    ++ic;
  }

  double t0=0,t1=0,t2=0,t3=0,t4=0,t5=0,t6=0;
  //
  int nx=n/npx;
  int ny=n/npy;
  int nz=(n/2+1)*2;
  int ny_half=ny/2;
  int ntot=nx*ny*nz*3;
  Timer out;
#pragma omp parallel reduction(+:t0,t1,t2,t3,t4,t5,t6)
  {
    container_type A(ntot),B(ntot);
    Vector<int> kx_map(nx);
    Vector<int> ky_map(ny_half);
    Vector<int> kz_map(nz);

    //this is arbitrary
    for(int i=0; i<kx_map.size(); ++i) kx_map[i]=(i<nx/2)?i:nx-i;
    for(int i=0; i<ky_map.size(); ++i) ky_map[i]=(i<ny/2)?i:ny-i;
    for(int i=0; i<kz_map.size(); ++i) kz_map[i]=(i<nz/2)?i:nz-i;

    real_type dt_xyz3=0.0, dt_3xyz=0.0, dt_xy3z=0.0, dt_cxyz3=0.0;
    Timer myclock;
    vector<double> dtl(8,0.0);
    for(int i=0; i<niters; ++i)
    {
      myclock.restart();
      dns_ops::vorticity_xyz3(A.data(),B.data()
          ,kx_map.data(),kx_map.size()
          ,ky_map.data(),ky_map.size()
          ,kz_map.data(),kz_map.size()
          );
      dtl[0]+=myclock.elapsed();

      int zero=0;
      int nxyz3=ntot/2;
#if defined(ENABLE_FORTRAN)
      myclock.restart();
      vxyz3(reinterpret_cast<complex_type*>(A.data())
          ,reinterpret_cast<complex_type*>(B.data())
          ,kx_map.data(),&nx
          ,ky_map.data(),&ny_half
          ,kz_map.data(),&nz
          ,&nxyz3);
      dtl[1]+=myclock.elapsed();
#endif

      myclock.restart();
      dns_ops::vorticity_xy3z(A.data(),B.data()
          ,kx_map.data(),kx_map.size()
          ,ky_map.data(),ky_map.size()
          ,kz_map.data(),kz_map.size()
          );
      dtl[2]+=myclock.elapsed();

#if defined(ENABLE_FORTRAN)
      myclock.restart();
      vxy3z(reinterpret_cast<complex_type*>(A.data())
          ,reinterpret_cast<complex_type*>(B.data())
          ,kx_map.data(),&nx
          ,ky_map.data(),&ny_half
          ,kz_map.data(),&nz
          );
      dtl[3]+=myclock.elapsed();
#endif


      myclock.restart();
      dns_ops::vorticity_3xyz(A.data(),B.data()
          ,kx_map.data(),kx_map.size()
          ,ky_map.data(),ky_map.size()
          ,kz_map.data(),kz_map.size()
          );
      dtl[4]+=myclock.elapsed();

#if defined(ENABLE_FORTRAN)
      myclock.restart();
      v3xyz(reinterpret_cast<complex_type*>(A.data())
          ,reinterpret_cast<complex_type*>(B.data())
          ,kx_map.data(),&nx
          ,ky_map.data(),&ny_half
          ,kz_map.data(),&nz
          );
      dtl[5]+=myclock.elapsed();
#endif


      myclock.restart();
      dns_ops::vorticity_xyz3(reinterpret_cast<complex_type*>(A.data())
          ,reinterpret_cast<complex_type*>(B.data())
          ,kx_map.data(),kx_map.size()
          ,ky_map.data(),ky_map.size()
          ,kz_map.data(),kz_map.size()
          );
      dtl[6]+=myclock.elapsed();
    }

    t0+=dtl[0];
    t1+=dtl[1];
    t2+=dtl[2];
    t3+=dtl[3];
    t4+=dtl[4];
    t5+=dtl[5];
    t6+=dtl[6];
  }

  {
    vector<double> timing(7,0.0);
    timing[0]=t0; timing[1]=t1; timing[2]=t2; timing[3]=t3;
    timing[0]=t4; timing[1]=t5; timing[2]=t6;

    mpi::reduce(mycomm,timing);

    real_type norm=1.0/static_cast<real_type>(niters*omp_get_max_threads()*mycomm.size());
    cout << "size npencil mpi openmp c-xyz3 f-xyz3 c-xy3z f-xy3z c-3xyz f-3xyz c-xyz3-complex\n";
    //cout << ntot << " " 
    cout << n << " "   << nx*ny << " "
      << mycomm.size() << " " << omp_get_max_threads() << " " 
      << t0*norm << " " << t1*norm << " " << t2*norm << " " 
      << t3*norm << " " << t4*norm << " " << t5*norm << " "  << t6*norm <<  " "
      << out.elapsed() << " vorticity " <<endl;
  }

  //ostringstream os;
  //os << "dims nz nfft tasks threads loops time bigtime " << endl;
  //os << n_var << " "
  //  << nz << " " << nc << " " 
  //  << mycomm.size() << " " << omp_get_max_threads() << " " << nb_tot << " "
  //  << dt_tot/static_cast<double>(np*niters) << " " 
  //  << node_clock.elapsed()/static_cast<double>(niters);
  //cout << os.str() << endl;
  return 0;
}
